clear all
clc
close all

load AALabData
Row=row;
Column=column;

load AANewData
Row=[Row;row'];
Column=[Column;column'];

load AAMTData
Row=[Row;row];
Column=[Column;column];

ROW=zeros(486,9);
for i=1:486
    ROW(i,:)=reshape(Row{i},[9,1])';
end
[~,I]=unique(ROW,'rows');

Row=Row(I,:);
Column=Column(I,:);

n=length(Row);

Q=100;
fakedistr=cell(1,Q);
ErrorUniform=zeros(1,Q);
ErrorLogitPCHM=zeros(1,Q);
complete=9999*ones(1,Q);
fittau=cell(1,Q);
fitlambda=cell(1,Q);
errparam=cell(1,Q);

testtau=0:0.01:0.1;
u=length(testtau);
lambda=0:0.01:0.1;
s=length(lambda);

indtau=zeros(1,Q);
indlambda=zeros(1,Q);

load FM

for q=1:Q
    
    % Generate fake distributions
    fakedistr{q}=genfake(Row,Dominant,Dominated);
    %fakedistr{q}=genfakeNoRestr(Row);
    Freq=fakedistr{q};
    
    % Predict uniform
    
    err=zeros(1,n);
    for i=1:n
        err(i)=Freq(i,1)*log((Freq(i,1)+eps)/(1/3))+Freq(i,2)*log((Freq(i,2)+eps)/(1/3))+Freq(i,3)*log((Freq(i,3)+eps)/(1/3));
    end
    ErrorUniform(q)=mean(err);
    
    % Run PCHM
    
    errparam{q}=zeros(u,s);
    for k=1:u
        for p=1:s
            errparam{q}(k,p)=CHfixedtaulambda(testtau(k),lambda(p),Row,Column,Freq);
        end
    end
    [min_val, position] = min(errparam{q}(:));
    [indtau(q),indlambda(q)]= ind2sub(size(errparam{q}),position);
    ErrorLogitPCHM(q)=min(min(errparam{q}));
 
    
end

figure
histogram(1-complete,4) 
set(gca,'fontsize', 20)
xlim([0.9, 1])
%xlim([min(1-complete),1])
saveas(gcf,'histogram-logitPCHM.png', 'png')

'Restrictiveness'
r=mean(ErrorLogitPCHM)/mean(ErrorUniform)

'Standard error'
sampleVarLPCHM = mean((ErrorLogitPCHM-mean(ErrorLogitPCHM)).^2);
sampleVarUniform = mean((ErrorUniform-mean(ErrorUniform)).^2);
covar = mean((ErrorLogitPCHM-mean(ErrorLogitPCHM)).*(ErrorUniform-mean(ErrorUniform)));

SE = sqrt((sampleVarLL1 - 2 * r * covar + r^2 * sampleVarUniform) / mean(ErrorUniform)^2) * (1/sqrt(Q))


